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Abstract 

In this correspondence the cumulants of the mutual information of the flat Rayleigh fading amplify-and-forward 
O ■ MIMO relay channel without direct link between source and destination are derived in the large array limit. The 

analysis is based on the replica trick and covers both spatially independent and correlated fading in the first and 
the second hop, while beamforming at all terminals is restricted to deterministic weight matrices. Expressions for 
mean and variance of the mutual information are obtained. Their parameters are determined by a nonlinear equation 
system. All higher cumulants are shown to vanish as the number of antennas n goes to infinity. In conclusion the 
distribution of the mutual information I becomes Gaussian in the large n limit and is completely characterized by 
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the expressions obtained for mean and variance of /. Comparisons with simulation results show that the asymptotic 



results serve as excellent approximations for systems with only few antennas at each node. The derivation of the 



results follows the technique formalized by Moustakas et al. in [1]. Although the evaluations are more involved for 
the MIMO relay channel compared to point-to-point MIMO channels, the structure of the results is surprisingly 
simple again. In particular an elegant formula for the mean of the mutual information is obtained, i.e., the ergodic 
capacity of the two-hop amplify-and-forward MIMO relay channel without direct link. 

Index Terms 

MIMO relay channel, amplify-and-forward, replica analysis, random matrix theory, large antenna number limit, 
cumulants of mutual information, correlated channels. 



I. Introduction 

Cooperative relaying has obtained major attention in the wireless communications community in recent 
years due to its various potentials regarding the enhancement of diversity, achievable rates and range. An 
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important milestone within the wide scope of this field is the understanding of the fundamental limits 
of the MIMO relay channel. Such a channel consists of a source, relay and destination terminal, each 
equipped with multiple antennas. 

Generally, there are different ways of including relays in the transmission between a source and a 
destination terminal. Most commonly relays are introduced to either decode the noisy signal from the 
source or another relay, to re-encoded the signal and to transmit it to another relay {multi-hop) or the 
destination terminal {two-hop). Or the relay simply forwards a linearly modified version of the noisy 
signal. These relaying strategies are referred to as decode-and forward (DF) and amplify-and-forward (AF), 
respectively. Currently, the simple AF approach seems to be promising in many practical applications, 
e.g., since it is power efficient, does not introduce any decoding delay and achieves full diversity. Another 
approach is the so called compress-and-forward strategy (CF), which quantizes the received signal and 
re-encodes the resulting samples efficiently. 

We briefly give an overview over important contributions to the field of cooperative communications and 
relaying. The capability of relays to provide diversity for combating multipath fading has been studied 
in [2], [3] and [4]. In [5] the potential of spatial multiplexing gain enhancement in correlated fading 
channels by means of relays has been demonstrated. Tight upper and lower bounds on the capacity of 
the fading relay channel are provided in [6], [7], [8], [9], and [10]. Furthermore, in [11] the capacity has 
been shown to scale like N log K for the fading MIMO relay channel , where N is the number of source 
and destination antennas and K is the number of relays. 

In this paper we focus on the two-hop amplify-and-forward MIMO relay channel with either i.i.d. or 
correlated Rayleigh fading channel matrices. Our quantities of interest are the cumulant moments of the 
mutual information of this channel. Of particular importance in this context are its mean and variance. 
While the mean completely determines the long term achievable rate in a fast fading communication 
channel, the variance is crucial for the characterization of the outage capacity of a channel, which is 
commonly the quantity of interest in slow fading channels. Seeking for closed form expressions of 
cumulant moments of the mutual information in MIMO systems usually is a hopeless task. For the 
conventional point-to-point MIMO channel it therefore turned out to be useful to defer the analysis to the 
regime of large antenna numbers. For the i.i.d. Rayleigh fading MIMO channel closed form expressions 
were obtained in [12] and [13]. For correlated fading at either transmitter or receiver side the mean was 
derived [14], and [15] finally provided the mean for the case of MIMO interference. All these results are 
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obtained via the deterministic asymptotic eigenvalue spectra of the respective matrices appearing in the 
capacity log det-formula. 

Higher moments were also considered, e.g., in [16], [17] and [1], where the distribution in the large 
antenna limit was identified to be Gaussian. Generally, these large array results turned out to be very 
tight approximations of the respective quantities in finite dimensional systems. For amplify-and-forward 
MIMO relay channels only little progress has been achieved so far even in the large array limit. The mean 
mutual information of Rayleigh fading amplify-and-forward MIMO relay channels in the large array limit 
has been studied in [18] for the special case of a forwarding matrix proportional to the identity matrix 
and uncorrelated channel matrices. In this paper a fourth order equation for the Stieltjes transform of the 
corresponding asymptotic eigenvalue spectrum is found, which allows for a numerical evaluation of the 
mean mutual information. Since even for this special case no analytic solution is possible, the classical 
approach of evaluating the mean mutual information via its asymptotic eigenvalue spectrum does not seem 
to be promising for the AF MIMO relay channels. 

The key tool enabling the evaluation of the cumulant moments of the mutual information in the large 
array limit in this paper is the so called replica method. It was introduced by Edwards and Anderson in 
[19] and has its origins in physics where it is applied to large random systems, as they arise, e.g., in 
statistical mechanics. In the context of channel capacity it was applied by Tanaka in [20] for the first time. 
Moustakas et al. [1] finally used a framework utilizing the replica trick developed in [21] to evaluate the 
cumulant moments of the mutual information of the Rayleigh fading MIMO channel in the presence of 
correlated interference. The paper [1] is formulated in a very explicatory way and this correspondence 
goes very much along the lines of this reference. Though not being proven in a rigorous way yet, the 
replica method is a particularly attractive tool when dealing with functions of large random matrices, 
since it allows for the evaluation of arbitrary moments. Free probability theory, e.g., only allows for the 
evaluation of the mean, e.g., [22]. There are also some large array results by Muller that are of importance 
for amplify-and-forward relay channels. He applied free probability theory to concatenated vector fading 
channels in [23] (two hops) and [24] (infinitely many hops), which can be considered as multi-hop MIMO 
channels with noiseless relays. The contributions of this paper are summarized as follows: 

• In the large array limit we derive mean and variance of the mutual information of the two-hop MIMO 
AF relay channel without direct link where the channel matrices are modelled as Kronecker correlated 
Rayleigh fading channels while the precoding matrix at the source and also the forwarding matrix at 
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the relay are deterministic and constant over time. The obtained expression depends on coefficients 
that are determined by a system of six nonlinear equations. 

• We show that all higher cumulant moments are C(n _1 ) or smaller and thus vanish as n grows large. 
Accordingly, we conclude that the mutual information is Gaussian distributed with mean and variance 
given by our derived expressions in the large n limit. 

• Considering that not all doubts about the replica method are dispelled yet, we verify the obtained 
expressions by means of computer simulations and thus confirm that the replica method indeed works 
out in our problem. 

II. The Channel and its Mutual Information 

The two-hop MIMO amplify-and-forward relay channel under consideration is defined as follows. Three 
terminals are equipped with n s (source), n T (relay), and n d (destination) antennas, respectively. We allow 
for communication from source to relay and from relay to destination. Particularly, we do not allot a direct 
communication link between source and destination. Both the uplink (first hop from source to relay) and 
the downlink (second hop from relay to destination) are modelled as frequency-flat, i.e., the transmit 
symbol duration is much longer than the delay spread of up- and downlink. We denote the channel matrix 
of the uplink by Hi G C nrXn % the one of the downlink by H 2 G C™ dXnr . Furthermore, we assume that 
the relays process the received signals linearly. The matrix performing this linear mapping is denoted 
F r G C nrXnr and called the "forwarding matrix" in the following. With s the transmit symbol vector, a 
precoding matrix F s G C nsXns and n r and n d the relay and destination noise vectors respectively, the 
end-to-end input-output-relation of this channel is then given by 

y = H 2 F r HiF s s + H 2 F r n r + n d . (1) 

The system is represented in a block diagram in Fig. \T\ 

The elements of the channel matrices Hi and H 2 will be assumed to be zero mean circular symmetric 
complex Gaussian (ZMCSCG) random variables with covariance matrices as defined in the Kronecker 
model [25]: 

E[vec(Hi)vec(Hi) H ] = Tj ® R r , such that Tr{T s } = n s and Tr{R r } = n r , (2) 
E[vec(H 2 )vec(H 2 ) H ] = ® R d , such that Tr{T s } = n r and Tr{R d } = n d , (3) 

where vec(X) stacks X into a vector columnwise, cg> denotes the Kronecker product while Tr(-) and 
(■) T denote the trace and transposition operator, respectively. T s G C nsXn % R r G C nrXnr , T r G C nrX ™ r 
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and R d G C ndXrid are the (positive definite) covariance matrices of the antenna arrays at the respective 
terminals. These matrices are required to have full rank for the analysis below. We remind the reader 
that matrices following Gaussian distributions defined by covariance matrices as in © and © can be 
generated from a spatially white matrix H„, - in our case through the mappings 

Hx = R?H W)1 Tf (4) 

and 

H 2 = R|H Mj1 TJ. (5) 

The above described correlation model thus uses separable correlations, which is a commonly accepted 
assumption for wireless MEMO channels. 

Since we will be confronted with products of covariance matrices later on, we need to introduce the 
operator (•)* for quadratic matrices, which zeroizes the smaller of two matrices in a product such that it 
adapts to the size of the other matrix, or leaves it untouched if the matrix is the bigger one. Thus, we 
ensure that products like R*R d are well defined. As long as two matrices A and B are Toeplitz-like a 
product A*B* always yields the same result irrespective of the corner(s) used for zeroising the smaller 
matrix. 

We assume all channel matrix elements to be constant during a certain interval and to change indepen- 
dently from interval to interval (block fading). The input symbols are chosen to be i.i.d. ZMCSCGs with 
variance p, i.e., E[ss H ] = p/n s l ns , the additive noise at relay and destination is assumed to be white in both 
space and time and is modelled as ZMCSCG with unit variance, i.e., E^nl 3 ] = I r and E[n d rij] = I nd . 

The assumptions on the channel state information (CSI) are as follows: The destination perfectly knows 
the instantaneous channel matrices Hi and H 2 as well as F s and F s . The source and the relay only know 
the second order statistics of Hi and H 2 , i.e., the corresponding covariance matrices. In particular this 
implies, that the forwarding matrix can only depend on the covariance matrices of Hi and H 2 , but not 
on the instantaneous channel realizations. Its elements thus are deterministic and remain constant over 
time. Our analysis could only capture time-varying forwarding matrices that are Gaussian. However, 
forwarding matrices chosen based on the current channel realization would not be Gaussian in general. It 
will be useful to decompose the forwarding matrix into a scaling factor ^Ja/n r and a matrix F r fulfilling 
Tr{F r Ff} = n r . We will denote a as the power gain of the forwarding matrix. 



tioifl 



With Tr{F s F^} = n s the mutual informationj conditioned on Hx and H 2 in nats per channel use can 
be written as 

'det (l nd + £H 2 F r F r H H« + ^H 2 F r HiF s FH H HFH H H 




det (l nd + ^H 2 F r FH H H) 



(6) 



I nd + — H 2 F r Ff (7) 
corresponds to the overall noise covariance matrix at destination and 

^H 2 F r HiF s FfH?FfH? (8) 

n s n r ' 

corresponds to the signal plus noise covariance matrix at the destination. Since the forwarding matrix 
does not depend on the instantaneous channel realizations by assumption, it can be incorporated into T r 
according to 

f r 4 F r T r Ff . (9) 
Similarly F s can be incorporated into T s according to 

f s 4 F s T s Ff . (10) 

Refer to the extended block diagram in Fig. [2] for an illustration. 

1 _ 1 _ 1 „ 1 

In terms of the respective equivalent channels Hx = R^H^xTs 2 and H 2 = R^H^^Tr © can be 

rewritten as 



det (l nd + ^H 2 H« + ^H 2 HxHH H ^ 

/(s;y) = ln ^ ^— '- . (11) 

det (id + ^H 2 HH 



In the subsequent sections we will work with (fTTj) and will drop the tildes again for the sake of clarity. 

Due to the randomness in Hx and H 2 also I is a random variable. The theorem stated in the following 
section will fully characterize the distribution of / in the limit of large antenna numbers. 

'in this chapter we pass on the common pre-log factor 1/2, which accounts for the use of two time slots necessary in half-duplex relay 
protocols. 



III. Results 

We formulate our results in the subsequent theorem. Whenever we use the notation 0(f(n)) in the 
following we assume that n s , n r and n& grow to infinity with all ratios among them fixed. 
Theorem 1: For the mutual information / as defined in (flTT) 

• the mean is 0(n) and given by 

E[J] = ln(det(I„ s +p Sl T s )) + ln(det(I nr + as 2 T r )) 

- In (det (I nr + as 3 T T )) - In (det (I„ d + t 3 R d )) 
+ In (det (I max(nr , nd) + t 2 R^ + M 2 R r *R*)) 

- (n s s 1 t 1 + n T s 2 t 2 - n T s 3 t 3 ) + 0{n~ l ) (12) 

with 

n s h = TrjpTs^+psxTs]- 1 } (13) 
nM = Ti{aT T [I ni + as 2 T I ]- 1 } (14) 
n v t 3 = TrjaTr^ + asaT^ 1 } (15) 
n sSl = Tr{t 2 R r *R* [I„ r + M 2 R r *R* + ^R*]" 1 } (16) 
n T s 2 = Tr{(R T + t 2 R*R* d )[l nt + t 1 t 2 R* T R* d + t 2 R* d }- 1 } (17) 
n T s 3 = Tr{R d [I„ r +t 3 Rd] _1 }, (18) 

• the variance is 0(1) and given by 

Var[/] = - ln(det(Vi)) - ln(det(V 2 )) + 2 ln(det(V 3 )) + 0{n" 2 ) (19) 





1 



























1 



















«? } 










„?> 


1 


„?> 






















1 












T,« 















«?> 


1 
























1 


u % 





v® 


u ll 




















1 










(*\ \ 

v { ' 
u 8 







1 M> 
u 9 







1 


vi 1] 


u ll 








„?> 

















1 


V 


4 







6 





11 





„s> 1 


7 / 








V 


2 = 


A, ( 2 ) 
I 1 


! ^ 


1 














( 1 




(3) 

-v[ 







(3) \ 
1 








v 3 = 







1 

(3) 
f [ 


-i 






(3) 












(3) 
\ V 2 





-v {3) - 


- 1 J 







^{(pT^I + p^T,]- 1 ) 2 } 
-^Tr{(«T r [I + a^Tr]- 1 ) 2 } 

-Tr |R r *R* ([I + t 2 R*] [I + t 2 R* d + ti^R^R*]- 1 ) 2 } 

-Tr {R r (t 2 R* [I + t 2 R* + M^R*]- 1 ) 2 } 

Tr { (t 2 R r *R^ [I + t 2 R^ + i^RM" 1 ) 2 } 

Tr { (*iR r *R^ [I + t 2 R* d ] [I + t 2 R^ + MaR^]" 1 ) 2 } 

Tr { (M 2 R r *Rf [I + t 2 R^ + M^R^ 1 ) 2 } 

Tr { ([R* + hhR^Rf ] [I + t 2 R* + M^R*]" 1 ) 2 } 

Tr |t 2 R*RjR^j [I + t 2 R*] ([I + t 2 R* + t^.R^R*]- 1 ) 2 } 

-Tr |t 2 R* [I + t 2 R*] (tiR r *R^ [I + t 2 R^ + ti^R^R^^ 1 ) 2 ] 
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= Tr|t!R r *R* 2 [I + t 2 R*] [I + tit 2 R r *Rd] 






x ([I + t^ + M^R*]- 1 ) | 


(33) 


u n 


= -Tr {M 2 R r *Rf [I + M 2 R r *R^] ([I + t 2 R^ + fctaRJRST 1 ) 2 } 


(34) 




1 f 1 

= ^{(aT^ + a^Tj- 1 ) 2 } 


(35) 




= Tr{(R d [I + t 3 Rd] _1 ) 2 } 


(36) 


(3) 

v{ 


i 

= -Tr {aR*T; [I + t 3 Rd] [I + a^TX 1 } 


(37) 


(3) 


= — Tr {tit 2 R^T*R*Rd [I + t 2 R* + tx^R^R*]" 1 [I + assT^ 1 } 


(38) 


(3) 


= — TrlR^T* [I + tit 2 R r *RS] 






x [I + t 2 R* + t^R^R*]" 1 [I + assT;}- 1 } 


(39) 


„(3) 


= — Tr{R*jT* [tiR* + tit 2 R*R d ] 






x [I + t 2 R* + t^R^R*]" 1 [I + assT^r 1 }, 


(40) 



• all higher cumulant moments (skewness, kurtosis, etc.) are 0(n _1 ) or smaller and thus vanish, as n 
grows large. Consequently, the mutual information / is Gaussian distributed random variable in the 
large n limit. 

IV. Mathematical Tools 

In this section we briefly repeat the mathematical tools we will use in the proof of the theorem. These 
are (cumulant) moment generating functions, the replica method and saddle point integration. At the same 
time we shall give a brief outline of the proof, which we will provide in full detail in Section |VJ 

A. Generating Functions 

We define the moment generating function of the mutual information / as follows: 

9l (u) = E [e-" 1 ] . (41) 

This definition differs from the standard definition in the sign of the argument of the exponential function. 
The minus sign used in the definition above will simplify our analysis later on. Provided that the moment 
generating function exists in an interval around i/ = 0we may expand (|4TI) into a series in the following 
way 

9l (u) = 1 - v • E[J] + % • E[/ 2 ] - £ ■ E[/ 3 ] + . . . . (42) 
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We will also consider the cumulant generating function of /, which is defined as ln<?j(V). Expanded into 
a Taylor series around zero it is given by 



v 2 .. „ AM' 



H 9l (u)) = -v ■ E[ln/(X)] + - ■ Var[ln/(X)] + £ ^^C p , (43) 

^ P =3 p - 

with C p the pth cumulant moment. Once we have found this series, it is thus easy to extract mean and 
variance by a simple comparison of coefficients. Furthermore, since a Gaussian random variable has the 
unique property that only a finite number of its cumulants are nonzero (more precisely its mean and 
variance), we will be able to proof the asymptotic Gaussianity of I by showing that the cumulants C p die 
out for p > 2 in the large n limit. 

B. Integral Identities 

We will need some useful integral identities in order to evaluate the moment generating function. Before 
stating them we introduce a compact notation for products of differentials arising when integration over 
elements of matrices is performed. With % = T as well as $IZ and ^sZ the real and imaginary part of 
a complex variable Z, we define the following integral measures, which are completely identical with the 
notation used in [1]: 

d c X = Y\ Y\. d^-Xijd^sXij, for X i3 complex variables, 



(44) 

d g X. = — Y\ Y\. dXijdX^ r, for Xjj and Grassmann variables, (45) 

i 3 

d/i(X, Y) = Y\ Y\. dXijdYji, for Xy and Yij complex variables. (46) 



2tt 

i j 



2m 

j 



The defining properties of a Grassmann variable are listed in Appendix IB With this notation as well as 
Cg> the Kronecker product operator we specify the following identities, which are all proven in [1]: 
. For M G C nxn , N e C vxu positive definite, and X, A, B e C nxv we have 

J exp (^-^Tr (NX H MX + A H X - X H B)^ d c X 

= (det (N ® M)) -1 exp (-^Tr (N -1 A H M~ 1 B) J . (47) 

• For M G C nxn , N G C uxu positive definite, and A, X and X, B n x v and uxn matrices, respectively, 
whose entries are Grassmann variables, we have 

J exp (Tr (NXMX + AX + XB)) rf 9 X 

= det (N <g> M) exp (Tr (N^AM^B) ) . (48) 
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. For X, Y, A, B e C uxu we have 

J exp (Tr (XY - XA - BY)) d/j(X, Y) = exp (-Tr (AB)) . (49) 

The application of these identities is known as the replica trick, which introduces multiple copies of the 
Gaussian integration that arises when computing the expectation of exp(— vT) over the elements of Hi 
and H 2 . We emphasize that the machinery of repeatedly applying the above identities in the evaluation 
of giiy) (see Section IV-AI ) requires v to be a positive integer. In order to extract the (cumulant) moments 
of I from the respective generating function we thus will need to assume that gi(u) can be analytically 
continued at least in the positive vicinity of zero in the end. This assumption will be applied without 
being proven anywhere in the literature yet. Nevertheless, all results obtained based on this assumption - 
including those derived below - show a perfect match with results obtained through computer simulations. 

C. Saddle Point Integration 

For the final evaluation of the moment generating function we will use the saddle point method. In its 
simplest form it is an useful tool to solve integrals of the form 

lim / e -«-*(*i.~.**) . dx x ■ ■ ■ dx k , (50) 

n— »oo / 

where \&(-, . . . , •) is some function with well defined Hessian at its global minimum. We will use it with 
a slightly different expression in this paper. For the sake of clarity we will consider the univariate case in 
this section. In the actual proof of the Theorem we will then deal with integrals over multiple variables. 
Suppose we can rewrite the moment generating function of I in the form (as done in Section IV-AI ) 

gi {v) = J e~ nx ' u ' n) dx. (51) 

If we expand /(•) into a Taylor series in x around its global minimum in x we can write 

gi{y) = e - -^ '^'™) • J e ~\f" { x o,y,n)(x-x ) 2 -\f'" {x ,v,n){x-x ) i + ... . ^ ^2) 

where the (■)' operator denotes derivation for x. From this expansion and our particular function /(•, u, n), 
which will be multivariate indeed, it will possible to show that (|52l ) evaluates to 

9l (p) = exp (-vn- ^(x ) + v 2 ■ £ 2 (x )) + ^ v k ■ O^- 1 ), (53) 

fc>3 

with and £ 2 ( - ) functions that we determine in Section IV-BI The fact that 

(54) 



= (-!)'• 



u=0 
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immediately reveals that the leading terms of mean and variance are determined by £1 and £ 2 > respectively. 
The O^n^ 1 ) scaling of the residual terms is proven in Section IV-CI Comparing higi(v) to the right hand 
side of (|43l ) will reveal the higher order cumulants to be at most 0(n~ l ). Remember that we obtained 
(|43l ) as a Taylor expansion around v — 0. We thus have implicitly assumed that the limit n — ► oo and 
v — ► can be interchanged. This assumption is noncritical and made without proof in this paper. 

In the subsequent sections we will apply this procedure in a multivariate framework. f(-,u,n) will then 
be a function of multiple matrices with a appropriately defined integration measures (cf. next subsection), 
which appear inside trace and determinant operators. We will make a symmetry assumption called the 
hypothesis of replica invariance, namely that all these matrices are proportional to the identity matrix at 
the global minimum of /(•, is, n). This assumption is justified in [26]. Therefore no proof is provided in 
this paper. 

We highlight that it is this saddle point method that makes the following derivations a large n approxi- 
mation. If we had another tool capable to solve the critical integral for finite n the whole procedure could 
also be applied to obtain nonasymptotic results. 

V. Proof 

The equations in the prooln are somewhat involved. In order to make the proof clearly laid out and 
more compact we therefore omit the channel covariance matrices and assume antenna arrays of size n at 
each terminal at first instance. Both covariance matrices and the possibly different antenna numbers can 
be easily reintroduced at the end of the proof. For the sake of clarity we structure the proof into three 
parts corresponding to the subsections below. 

A. Applying the Replica Trick 

We introduce the auxiliary variables X, Y, Z,Wi,W 2 G C nxu and A, B, AB (u x n and n x v 
Grassmann matrices) and evaluate the moment generating function of / by means of identities (|47l) - ( |49l ) 

2 This proof is only rigorous in the case that the analytic continuation of giiy) to zero is indeed possible. Proving this in turn is a current 
research topic in mathematics. 
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x exp ( --Tr ( -R 5 X H X - Q 5 W"Z + -R 6 X H X + Q 6 Z H Z )) 

\ 2 \n n / J 

x exp ( --Tr f--BZQ 4 - 2Q 4 X H A + -AXR 7 + 2R 7 W^b') 

\ 2 \ n n J 

x exp ( --Tr f--AXQ 7 - 2Q 7 Z H B - 2-R 8 BB - 2Q 8 AA )) 

\ 2 \ n n J J 



J exp(— S) ■ dX, 



xd c W 1 d c W 2 d c Xd c Yd c Zd g Ad g BdX (59) 

(60) 



where we have combined the integral measures over the various Rj's an Q/s into the single integral 
measure 

dX 4 d / i(R 1 ,Q 1 )rf/i(R 2 ,Q 2 )rf/i(R3,Q 3 )rf/i(R 5 ,Q5) 

xcfyt(R 6 , Q 6 )d/i(R 8 , Q 8 ) ■ dR 4 dQ4dK 7 dQ 7 . (61) 

In (|56l ) we have firstly applied (|47l ) and (|48l ) (backwards) with M the argument of the determinant in 
nominator and denominator respectively, N = l uxu and A = B = nx „ in order to get rid of the 
determinants. Afterwards we again applied (l47l) (backwards) twice with M = I nxn and N = l uxu in 
order to split the products H 2 HiH^HJ? and H 2 HJ? at the expense of the introduced auxiliary matrices 
Y and Z. For the first application we have A = B = H^HJ?, for the second one A = B = H^. Exactly 
the same is done in (1571) again where we also break up the products H 2 Hi and H^H^ 1 . In (l58l) we get 
rid of the integrals over Hi and H 2 by twice applying (l47l) (forwards). In (l59l we split all quartic terms 
into quadratic terms by making use of (l48l) and (|49l) . We can get rid of all integrals but the outer one, by 
(forwards) applying identities (PTTT) and (l48l) again, and after some algebraic effort we obtain S as 



S = -Tr (RiQi + R 2 Q 2 + R 3 Q 3 + R4R4 + R5Q5 
+R 6 Q 6 — Q4Q4 + R7R7 — Q7Q7 — RsQs) 

+n In det (l„ + -Ri) + n In det (l u + -(R 2 + R 3 + R 5 + R 6 

V n / V n 

+n In det (I v + QiQ 2 ) + n In det (l„ - (I„ + QiQ 2 ) _1 Q5Q1Q3 + Qe) 

-n In det (l v + -Rs - - (I„ + Q1Q2) R7Q1R4 

V n n 

+- [I v - (I, + QiQ 2 ) 1 Q5Q1Q3 + Qe] " 



Oi r T / T ^ ^ 1 ^ ^ ^ ^ n-l 

n 

x [-R 7 QiQ 3 (I, + QiQ 2 ) _1 Q5Q1R4 (lu + Q1Q2)" 1 



+R 7 Q!Q3Q 4 (lu + Q1Q2) 1 + Q 7 (lu + Q1Q2) 1 Q5Q1R4 - Q 7 Q 4 ] ) 
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a 1 - 1 

I, + -(R 2 + R 3 + R 5 + R 6 ) 

n 

x [-R 4 R 7 + Q 4 R 7 + R 4 Q 7 - Q 4 Q 7 ] ^ . (62) 

At this point we have shaped the problem into the form of (l52l) . where the role of x is played by the 
introduced u x v auxiliary matrices. Note that there appears no matrix with one of its dimension equal 
to n in S anymore. 



-n In det ( \ v + Q 



a 

+ - 

n 



B. Evaluating Mean and Variance 

In order to evaluate the last remaining integral in (l60l) by means of saddle point integration we need to 
expand S into a Taylor series in <5Ri, <5Qi, . . . , <5R 8 , SQ 8 around its minimum. This expansion corresponds 
to the expansion in x in Section HV-Cl With S p denoting the pth order term in the series the expansion 
looks as follows 

S = S + S 2 + S 3 + ... (63) 

By symmetry all complex matrices are assumed to be proportional to the identity matrix at the minimum 
of S (replica symmetry), the Grassmann matrices have to vanish in order to obtain a real solution (by 
definition real numbers cannot be Grassmann numbers, since they commute). Thus, to develop the Taylor 
series (|63l) in this point we write 



Ri 


= r\nl v 


+ <5Ri 


(64) 


Qi 


= q x l v + £Qi 


(74) 


R 2 


= r 2 nl u 


+ £R 2 


(65) 


Q 2 


= q 2 l v + SQ 2 


(75) 


R 3 


= r z n\ v 


+ <5R 3 


(66) 


Q 3 


= q z \ v + ^Q 3 


(76) 


R 4 


= <5R 4 




(67) 


Q 4 


= (5Q 4 


(77) 


R 4 


= <5R 4 




(68) 


Q 4 


= 5Q A 


(78) 


R 5 


= r 5 nl v 


+ <5R 5 


(69) 


Q 5 


= q 5 l u + SQ 5 


(79) 


R 6 


= r§n\ v 


+ SR G 


(70) 


Qe 


= q % \ v + <5Q 6 


(80) 


Rr 


= <5R 7 




(71) 


Q 7 


= SQ 7 


(81) 


R 7 


= SK 7 




(72) 


Q 7 


= 6Q 7 


(82) 


Rs 


= r 8 nl u 


+ 6R S 


(73) 


Qs 


= q&l u + <5Q 8 


(83) 
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By definition S is given by (1621 evaluated at the minimum of S, i.e., 

S Q = v ■ n ■ {In (1 + pri) + In (1 + a(r 2 + r 3 + r 5 + r 6 )) 

+ In (1 + gig 2 - <?ig 3 <?5<?2 + q 6 + 9i929e) - In (1 + ar 8 ) 

- In (1 + g 8 ) - (ngi + r 2 q 2 + r 3 q 3 + r 5 g 5 + r 6 g 6 - r 8 g 8 )}. (84) 



The respective coefficients and have to ensure that <Si = 0. They are found by differentiating (1841) 
for each of them and setting the resulting expressions to zero. The derivatives for the r/s (note that we 
can summarize r 2 + r 3 + r 5 + r 6 + r$ = f 2 by symmetry) yield 

= qi ~~~ (85) = q 5 ~—^-r- (88) 

1 + pri 1 + ar 2 

= <?2-— — (86) = q e ~— — (89) 

1 + ar 2 1 + ar 2 

= q 3 -—^- (87) = g8 -_^_. (90) 

1 + ar 2 1 + ar s 

We see that q 2 = q 3 = q 5 = q 6 . Taking this into account the derivatives for the g/s yield 

o = n- q * (91) 
1 + gi?2 + <?2 

= r 2 - — ■ (92) 

1 + <?1<?2 + <?2 

= r 8 --^— . (93) 

The leading term thus simplifies to 

S = v ■ n ■ {In (1 + pri) + In (1 + ar 2 ) + In (1 + q 2 + 91^2) 

- In (1 + or 8 ) - In (1 + q 8 ) - (ngi + f 2 g 2 - r 8 g 8 )} = 1/ • ra • £1 (94) 

with 



* = rr^ <95) ri = + * <98) 

« = i - (96) f 2 = 1 + » (99) 

1 + ar 2 1 + gi?2 + <?2 

a 1 

58 = 7- (97) r = . (100) 

We note that £i(-) in (l94b is the multivariate version of the function mentioned in Section HV-C1 We see 

that £i(-) is 0(n°) and thus n • which will turn out to correspond to the mean of I in the large n 
limit, is 0(n). 
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At this point, we make use of the variable transformations H x — ► 6H X and Q x — > <5Q a . for x = 1...8, 
which preserve the integral measures. We indicate that transformation by denoting the respective integral 
measure. Furthermore, we define 



X ab ~ [^-^l,o6, 8R2,ab, SR^ a b, SR^ab, ^Re,ab, <%?l,a&, <5Q2,aft, $Q?,,ab-, ^Q^,ab-, 8Qe,< 
X aA ~ [^Rs,ab,SQ8,ab} T 

X a6 ~ [^4,a&, SR^ab-, SQ^ab-, ab , SRj^b, SRr,ab, <5<?7,a6, a fo] T - 



(101) 
(102) 
(103) 



With this notation we can write the moment generating function in terms of the Hessians of (1621 ), V 1; V 2 
and V 3 , as defined in d20j) - d22l) as 



9i\y) 



-So 



-So 



...)-d\ 



J exp(-S 2 - S 3 ~ S 4 
■ J exp(-S 2 ) • jl - [S 3 + S 4 + ...} + i [S 3 + S 4 + . . .] 

5EE x 2 Tv ' x i 



(104) 
dA (105) 



i=l a, 6=1 



x ^ 1 - [S 3 + <S 4 + . . .] + i [S 3 + S 4 + • • -f 



(106) 



-So 



det Vi det V 2 



(detV 3 ) 



/-pi~EE x f v ^ 



i=l a, 6=1 



•[S3 + S4 



2[S 3 + S 4 



■dX 



(107) 



In (11051 ) we expanded exp(— S 3 — S 4 — . . .) into a series. The evaluation of the integral over the first term 
in (11061 ) is provided in [1]. We note that £ 2 (-) — - In det |Vi| - In | det V 2 | + 2 In det |V 3 |, which will 
turn out to correspond to the variance of / in the large n limit, is 0(1). Again, £ 2 (-) is the multivariate 
version of the function mentioned in Section HV-Cl 



C. Proving Gaussianity 

We will next show, that the remaining integral expression 

3 v 



i=l a,b=l 



[S 3 + S 4 + ...] + - [S 3 + S 4 + . . .] 2 



■ rfA (108) 



is 0(n x ). To see this we need to consider the various Taylor coefficients of the S p for p > 2 first. By 
inspecting (|62|) we note that 
la) a differentiation for either R ls R 2 , R 3 , R 5 , R 6 or R 8 yields a multiplication by a factor 1/n, 
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2a) a differentiation for either Q l5 Q 2 , Q3, Q5, Q6 or Q 8 does not change the order with respect to n, 
3a) two differentiations for Grassmann variables (note that odd numbers of differentiations yield zero 
Taylor coefficients) yield a multiplication by a factor \/n. 
Accordingly a Taylor coefficient resulting from i, j and k differentiations of the first, second and third 
type, respectively will be 0(n 1 ~ t ~ k / 2 ). Also, a product of t Taylor coefficients resulting from i ± , j\, hi, 
32, k 2 , . . ., i t , it, h differentiations of the first, second and third type, each, will be C^n* - ^^^'/ 2 )). 
Next, consider integrals of the form 

/ ex PhEE x 2 Tv ' x 2 II II SC *r II SR kl SQ k2 SQ k3 SR k4 -d\. (109) 

For the complex matrices Wick's theorem allows us to split the integral into sums of products of integrals 
involving only quadratic correlations. Furthermore, it states that for odd numbers of factors the integral 
evaluates to zero. Ignoring the Grassmann matrices for the moment we can extract the order of these 
correlations in the following. We define V as the joint Hessian 



V=| I (110) 

v 2 

and note that det(V) is 0(1). Also, we define x = [x^' , x^ 6 ; ' ] T and denote the integral measure d\ 
without all Grassmann contributions by d\. With this notation we can extract the orders of the three kinds 
of arising quadratic correlations by applying the second part of Wick's theorem: 
lb) 



Vi 



/ ex P ( ^2 xTVx ) ' $ R i,ab ■ SRj,cd ■ d\ 

J V 0,6=1 / 



= « c |det(V)|-- • l det(v) >- = 0(n), (111) 



2b) 

J eX P f ~\ ^2 xTVx J ■ 5 Qi,ab ■ $Qj,cd ■ d\ 



a, 6=1 



SaM det(V)|— • Qet d l J (v) } = 0(n- 1 ), (112) 

(113) 
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3b) 



exp 



1 v 

In 



x T Vx • 5R i)CLb ■ 5Qj, cd ■ dX 



o,b=l 



-5 Ac|det(V)| 



det(V( 2i " 1 



0(1). 



(114) 



det(V) 

By det(V^ a ' b ^) we denote the sub-determinant when the ath row and the 6th column in the matrix is 
deleted, 5 xy denotes the Kronecker delta function. The orders follow, since deleting odd lines/columns in 
V amounts to a multiplication of the respective determinant by a factor which is 0(n), while deleting 
even lines/columns in V amounts to a multiplication of the respective determinant by a factor which is 
0(n~ l ). The Grassmannian integrations are easily verified to yield 0(n°) factors, since also the elements 
of V 3 are O(n ). 

Combining la) and lb), 2a) and 2b) as well as 3a) and 3b), we can finally summarize, that terms 
resulting from the evaluation of (11081) are 



„ / ( ^' 'x+jx+kx \ V— V . ... 

O [n ^«=i 2 ) , if y J x + j x + k x is even 



x=l 



or zero otherwise. Here, t denotes the number of involved Taylor coefficients, i x ,j x ,k x the number of 
derivations of kind 1, 2 and 3. Note that i x ,j x and k x also correspond to the number of factors arising 



with the Taylor coefficient in the correlation. Since XL=i 
appearing terms in the integral are O^n^ 1 ) or smaller. 
We can thus rewrite (11071) as 

]detVidetV 2 



> t for p > 2, we conclude that all 



9i\y 



(detV 3 ) 



(115) 



After factoring out the determinant the cumulant generating function is given by 



In g^v) 



In < e 



-s ti 



det Vi det V 2 



(detV 3 ) 2 



det Vi det V 2 



(detV 3 ) 2 




(116) 



-v-n-b- — (In | det(V x )| + In | det(V 2 )| 

-2 In | det(V 3 )|) + ln(l + O^ 1 )) 



v 



-v ■ n ■ £i + "— ■ 6 + Oin- 1 ) 



(117) 
(118) 



A coefficient comparison with (1431) immediately reveals 

E[I]=n^ + 0(n- 1 ), 



(119) 
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and 

Var[J] =& + 0(n- 1 ). (120) 

Also the C p for p > 2 are and thus vanish for n — > oo. This implies that / is Gaussian distributed 

in this limit. Note, that indeed the residual term of the variance can be shown to be 0(n~ 2 ) in the same 
way as it is done in [1]. The reason behind this is that no 0(n _1 ) term proportional to u 2 is generated 
in (11071) . We skip this (in the present case very tedious) derivation for reasons of brevity. 

D. Reintroducing Covariance Matrices 

Finally, we reintroduce the omitted covariance matrices T s , R r , T r , Rd. In (1571) we see that the 

i i 

covariance matrices could be attached to the introduced auxiliary matrices as follows: Y H T S 2 , R r 2 W 2 , 

W^Tf, R|X, X H R|, TlWi ,WfR r , TjY, Z H T|, RfX, X H R|, TjZ, BT?, R|A, AR| and 

i 

T r 2 B. In (|59| ), we always obtain products involving only identical (square roots of) covariance matrices 
as factors. Thus, we can attach a factor T s to R l5 a factor R r to Q l5 factors T r to Q 2 , Q3, Q5, Q6, Rs> 
R 4 , Q 4 , R 7 and Q 7 , and factors R d to R 2 , R3, R5, R6, Qs> R4, Qa, and Q 7 . In (1621) these factors 
are combined in outer products, while the factor of n is removed and the \ v are replaced by I u . n . It is then 
obvious, that (|94|) translates to (fl~2)) . and also the entries of the Hessians (|20l) - ([22]) follow immediately. 
From the dimension of the covariance matrices we can now also conclude the respective antenna array 
dimension and thus also replace the n by either n d , n v or n d again. 

VI. Comparison with Simulation Results 

We verify the results stated in the theorem by means of computer experiments. For the mean this is done 
through Monte Carlo simulations. The respective plot is shown in Fig. |3j where we present the ergodic 
mutual information versus the SNR for n = n s = n T = n d = 2, 4 and 8. We observe that even for only two 
antennas the approximation is reasonable, for four antennas the match is close to perfect, while for eight 
antennas no difference between analytic approximation and numeric evaluation can be seen anymore. In 
order to also verify our results for the higher cumulant moments we compare the empirical cumulative 
distribution function (CDF) of the mutual information to a Gaussian CDF with mean and variance given 
in the theorem. The respective plot is shown in Fig. HI Again, we observe that the analytic approximation 
becomes tight indeed as n = n s = n r = n d increases. For n = 8 even the tails of the distribution are 
reasonably approximated, which is an important issue for the characterization of the outage capacity. Our 
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simulation results thus also demonstrate that the replica method - despite its deficiency of not being 
mathematically rigorous yet - indeed reveals the correct solution to our problem. 

VII. Conclusion 

Using the framework developped in [21] and [1] we evaluated the cumulant moments of the mutual 
information for MIMO amplify and forward relay channels in the asymptotic regime of large antenna 
numbers. Similarly to the case of ordinary point-to-point MIMO channels, we observe that all cumulant 
moments of order larger than two vanish as the antenna array sizes grow large and conclude that the 
respective mutual information is Gaussian distributed. For mean and variance we obtain expressions 
that allow for an analytic evaluation. Computer experiments show, that the derived expressions serve as 
excellent approximations even for channels with only very few antennas. The results confirm the linear 
scaling of the ergodic mutual information (0(n)) in the antenna array size and also reveal that the 
respective variance is 0{\) in the antenna number. 
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Appendix I 
Preliminaries of Grassmann Variables 

Grassmann algebra is a concept from mathematical physics. A Grassmann variable (also called an 
anticommuting number) is a quantity that anticommutes with other Grassmann numbers but commutes 
with (ordinary) complex numbers. With 0i, 9 2 Grassmann variables and A a complex number the defining 
properties are 

A0i = 0iA (121) 
9 1 9 2 = -0 2 0x . (122) 

With 3 another Grassmann variable further properties are 

0i(0 2 3 ) = 3 (0i0 2 ) (123) 
9\ = (124) 
exp(0 1 2 ) = l + 0!0 2 . (125) 
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Integration over Grassmann variables is defined by the following to properties 

dO = (126) 
6d6 = 1. (127) 



Note that also the differentials are anticommuting, i.e., d9id9 2 = —dBidB\, Further details about integrals 
over Grassmann variables such as variable transformation can be found in the Appendix of [1]. 

Appendix II 
Wick's Theorem 

With V G C NxN , x G C Nxl and an integral measure da(x) = l/y/2ndx\, . . . , dxN we have 

(det V)a J exp ^-^x T Vx^ ]Jx k da(x) (128) 

= ^(det VO' / exp f-I x T Vx J • • ar i(2 • da(x) • (129) 

pairs ' ' 

x... (130) 
x(det V)^ J exp ^x T Vxj ■ x iM _ 1 • x iM ■ dof(x) (131) 

if M is even. For odd M the expression evaluates to zero. The sum in (11281) is over all possible 
rearrangements of the orderings of the indexes such that different indexes are paired with each other 
(with each distinct pairing being counted once). 
Furthermore, we have that 

(det V)^ j exp ^x T Vxj • x iXj ■ da(x) = [V' 1 )^, (132) 

with [V -1 ]^ the element in the ith row and jth column of V~ x . We will also need that 

[V-\j = det V^'V det V (133) 

with V^ ,J ' an N — 1 x N — 1 matrix, where the ith row and the jth column of V^ 1 are deleted. 
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Fig. 1. Block diagram of the channel. 
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Fig. 2. Extended block diagram of the channel. 
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Fig. 3. Mutual information versus SNR for n = n s = n r = and i.i.d. channel matrix entries - solid lines are analytical approximations, 
circles, squares and diamonds mark true mutual information as obtained through Monte Carlo simulations. 
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Fig. 4. Cumulative distribution function of mutual information for n = n a = n r = rid and i.i.d. channel matrix entries. Dashed lines 
represent Gaussian distributions with analytically computed mean and variance. The solid lines are the empirical distributions obtained 
through simulations. 



